library(plotly)
library(reshape2)

set.seed(4500)

# x <- rnorm(50, mean = 5, sd = 0.7)
# 
# mu_vals <- seq(0, 10, length.out = 100)
# sigma_vals <- seq(0.2, 1.2, length.out = 100)
# 
# x <- rnorm(40, mean = 5, sd = 2.2)

# prod(dnorm(obv, mu, sd))


# set.seed(123)
# x <- rnorm(40, mean = 5, sd = 2.2)

x <- c(3.6, 4.7, 5.8)

# 2. 构建参数网格
mu_vals <- seq(2, 8, length.out = 300)
sigma_vals <- seq(0, 4, length.out = 300)

log_likelihood <- matrix(NA, nrow = length(mu_vals), ncol = length(sigma_vals))

n <- length(x)

for (i in seq_along(mu_vals)) {
  for (j in seq_along(sigma_vals)) {
    mu <- mu_vals[i]
    sigma <- sigma_vals[j]
    log_likelihood[i, j] <- exp(-n * log(sigma) - n/2 * log(2*pi) - sum((x - mu)^2) / (2 * sigma^2))
    # log_likelihood[i, j] <- prod(dnorm(x, mean = mu, sd = sigma))
  }
}


plot_ly(x = mu_vals,
        y = sigma_vals,
        z = log_likelihood,
        type = 'surface') %>%
  layout(title = "Normal log-likelihood Surface",
         scene = list(
           xaxis = list(title = 'μ'),
           yaxis = list(title = 'σ'),
           zaxis = list(title = 'log-likelihood')
         ))
# z <- log_likelihood
#
# dimnames(z) <- list(mu_vals, sigma_vals)
#
# df <- melt(z, varnames = c("mu", "sigma"), value.name = 'logLik')
#
# plot_ly(df,
#         x =  ~ mu,
#         y =  ~ sigma,
#         z =  ~ logLik) %>%
#   add_surface() %>%
#   layout(title = "Normal log-likelihood Surface",
#          scene = list(
#            xaxis = list(title = 'μ'),
#            yaxis = list(title = 'σ'),
#            zaxis = list(title = 'log-likelihood')
#          ))
